#!/bin/bash
set -e

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------




# data_thread_num=10
# cpu_thread_num=4

# java_memory=-Xmx16g
# java_memory=16g
echo java memory: $java_memory 

# snpeff_db_version=GRCh37.75 
echo $snpeff_db_version 

# header=\
# @HD\tVN:1.4\tGO:none\tSO:coordinate

vcf_path=${data_path}/vcf/gatk/

echo ref genome is: $ref_genome

# java -jar snpEff.jar download -v GRCh37.71


# 经 gatk 和 snpeff 处理后，anno FILTER 报错，所以先运行
# ${script_path}anno_fp_0424.sh 1.snp.vcf 1.snp.fp.vcf 
# ${script_path}anno_fp_0424.sh 1.indel.vcf 1.indel.fp.vcf 

echo 
echo 
echo gatk VariantFiltration 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.snp.filtered.vcf \
	--variant 1.snp.vcf \
	--filterExpression "FS > 60.0" \
	--filterName "FS_more_than_60" \
	--filterExpression "MQ < 40.0" \
	--filterName "MQ_less_than_40" \
	--filterExpression "HaplotypeScore > 200.0" \
	--filterName "HaplotypeScore_more_than_200"
	
echo 
echo 
echo gatk VariantFiltration 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.indel.filtered.vcf \
	--variant 1.indel.vcf \
	--filterExpression "FS > 200.0" \
	--filterName "FS_more_than_200"
	

echo 
echo 
echo gatk CombineVariants 
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T CombineVariants \
	--variant 1.snp.filtered.vcf \
	--variant 1.indel.filtered.vcf \
	-o 1.snp.indel.filtered.vcf \
	-genotypeMergeOptions UNSORTED
	


# Annotate using ENCODE's data:
# java -Xmx4g -jar snpEff.jar -v -interval db/encode/wgEncodeDukeDnase8988T.fdr01peaks.hg19.bb GRCh37.75 examples/test.1KG.vcf > test.1KG.ann_encode.vcf


echo 
echo 
echo snpsift annotate dbsnp
java -Xmx$java_memory -jar $snpsift \
annotate \
$data_path/ncbi/dbsnp/All_20150605.vcf.gz \
1.snp.indel.filtered.vcf \
> 1.snpsift.dbsnp.vcf

echo 
echo 
echo snpsift annotate dbsnp papu
java -Xmx$java_memory -jar $snpsift \
annotate \
$data_path/ncbi/dbsnp/All_20150605_papu.vcf.gz \
1.snpsift.dbsnp.vcf \
> 1.dbsnp.papu.vcf

echo 
echo 
echo snpsift annotate clinvar
java -Xmx$java_memory -jar $snpsift \
annotate \
$data_path/ncbi/clinvar/clinvar_20150901.vcf.gz \
1.dbsnp.papu.vcf \
> 1.dbsnp.papu.clinvar.vcf

echo 
echo 
echo snpsift annotate clinvar papu
java -Xmx$java_memory -jar $snpsift \
annotate \
$data_path/ncbi/clinvar/clinvar_20150901_papu.vcf.gz \
1.dbsnp.papu.clinvar.vcf \
> 1.clinvar.papu.vcf


echo 
echo 
echo snpsift annotate cosmic 
java -Xmx$java_memory -jar $snpsift \
annotate \
${data_path}/cosmic/CosmicCodingMuts.vcf \
1.clinvar.papu.vcf \
> 1.clinvar.cosmic.vcf


echo 
echo 
echo snpsift annotate varType 
java -Xmx$java_memory -jar $snpsift \
varType \
1.clinvar.cosmic.vcf \
> 1.cosmic.var_type.vcf


echo 
echo 
echo snpsift annotate gwasCat 
java -Xmx$java_memory -jar $snpsift \
gwasCat \
-db ${data_path}/snpeff/gwascatalog.txt \
1.cosmic.var_type.vcf \
> 1.var_type.gwas.vcf


echo 
echo 
echo snpsift dbnsfp 
java -Xmx$java_memory -jar $snpsift \
dbnsfp \
-db ${data_path}/dbnsfp/dbNSFP2.9.txt.gz \
-f Ensembl_transcriptid,Uniprot_acc,Interpro_domain,SIFT_score,Polyphen2_HVAR_score \
1.var_type.gwas.vcf \
> 1.gwas.dbnsfp.vcf


# echo 
# echo 
# echo $snpeff ann snpeff2gatk 
# time \
# java -Xmx$java_memory -jar $snpeff \
# -c ${snpeff_path}/snpEff.config \
# -v \
# -lof \
# -noStats \
# -o gatk \
# $snpeff_db_version \
# 1.gwas.dbnsfp.vcf \
# > 1.dbnsfp.snpeff2gatk.vcf \
# 
# -t 8 \ ??
	
	
	
echo 
echo 
echo snpeff ann snpeff2gatk 
java -Xmx$java_memory -jar $snpeff \
-c ${snpeff_path}/snpEff.config \
-v \
-lof \
-noStats \
$snpeff_db_version \
1.gwas.dbnsfp.vcf \
> 1.dbnsfp.snpeff2gatk.vcf
	


:<<!
echo
echo
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.gatk.filter.vcf \
	--variant 1.dbnsfp.snpeff2gatk.vcf \
	--filterExpression "QD < 2.0" \
	--filterName "gatk_QD_less_than_2" \
	--filterExpression "MQRankSum < -12.5" \
	--filterName "gatk_MQRankSum_less_than_-12.5" \
	--filterExpression "ReadPosRankSum < -8.0" \
	--filterName "gatk_ReadPosRankSum_less_than_-8" \
	--filterExpression "ReadPosRankSum < -20.0" \
	--filterName "gatk_indel_ReadPosRankSum_less_than_-20" \
	--filterExpression "InbreedingCoeff < -0.8" \
	--filterName "gatk_indel_InbreedingCoeff_less_than_-0.8"
	
!



cat 1.dbnsfp.snpeff2gatk.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $1","$2","$3","$4}else{print $0}}' \
> 1.header_modified.vcf

echo 
echo
echo snpsift filter 
java -Xmx$java_memory -jar $snpsift \
filter \
"(FILTER = 'PASS')" \
1.header_modified.vcf \
> 1.high0.vcf



echo
echo
echo gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
-R $ref_genome \
-T VariantFiltration \
-o 1.high.vcf \
--variant 1.high0.vcf \
--filterExpression "QD < 2.0" \
--filterName "gatk_QD_less_than_2" \
--filterExpression "MQRankSum < -12.5" \
--filterName "gatk_MQRankSum_less_than_-12.5" \
--filterExpression "ReadPosRankSum < -8.0" \
--filterName "gatk_ReadPosRankSum_less_than_-8" \
--filterExpression "ReadPosRankSum < -20.0" \
--filterName "gatk_indel_ReadPosRankSum_less_than_-20"

# --filterExpression "InbreedingCoeff < -0.8" \
# --filterName "gatk_indel_InbreedingCoeff_less_than_-0.8"



# echo 
# echo 
# echo $snpsift filter 
# time \
# java -Xmx$java_memory -jar $snpsift \
# filter \
# "(FILTER != 'PASS')" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 



$bgzip -c 1.header_modified.vcf > 1.header_modified.vcf.gz
$tabix 1.header_modified.vcf.gz

echo 
echo 
echo bcftools view for filter FILTER 
$bcftools view \
-f "FP_VARIANTCALLER,FP_PRIMERDIMER" \
1.header_modified.vcf.gz \
| bgzip -c \
> 1.fp.vcf.gz


cat 1.header_modified.vcf |awk '{if(/^#/){print}else{if(/[(^PASS|^FP)]/){print}}}' > 1.low.vcf


# gunzip -c 1.fp.vcf.gz > 1.fp.vcf
# $tabix 1.fp.vcf.gz

# $bgzip -c 1.high.vcf > 1.high.vcf.gz
# $tabix 1.high.vcf.gz
# $bcftools concat -aD 1.high.vcf.gz 1.fp.vcf.gz |bgzip -c > 1.high_fp.vcf.gz
# $tabix 1.high_fp.vcf.gz

# if [ ! -d isec ]; then
	# mkdir isec
# fi
# $bcftools isec -p ./isec 1.header_modified.vcf.gz 1.high_fp.vcf.gz

# cp ./isec/0000.vcf 1.low.vcf



# echo 
# echo 
# echo bcftools view for filter FILTER 
# time \
# $bcftools view \
# -f "FS_more_than_60,MQ_less_than_40,HaplotypeScore_more_than_200,FS_more_than_200" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 



# java -jar SnpSift.jar rmInfo test.snpeff.vcf EFF
